extractorRData <- function(file, object) {
#' Function for extracting an object from a .RData file created by R's save() command
#' Inputs: RData file, object name
E <- new.env()
load(file=file, envir=E)
return(get(object, envir=E, inherits=F))
}
# Load genotypes and phenotypes
load(paste0(data.dir,"ROSMAP_SV_4_association.RData"))
pheno_list_files = list.files(paste0(data.dir,"/",studyDataID,"/association_results"), pattern = "*.assoc.RData", full.names = T)
pheno_list = pheno_list_files[gsub("(.*)/(.*)\\.assoc\\.RData","\\2",pheno_list_files) %in% names(pheno_list)]
res.all = data.frame()
for (pheno_i in pheno_list){
#pheno_i = pheno_list[1]
print(pheno_i)
res = extractorRData(pheno_i, "assoc")
res_i = res[,c("rs.id","beta","SE", "pval")]
res_i$bonf = p.adjust(res_i$pval, method = "bonferroni")
res_i$fdr = p.adjust(res_i$pval, method = "fdr")
colnames(res_i) = c("sv_id","Estimate","Std. Error", "nom_p", "bonf" , "fdr")
res_i$pheno = gsub("(.*)\\/(.*)\\.assoc\\.RData","\\2",pheno_i)
res.all <- bind_rows(res.all,res_i)
}
vcf.meta$sv_info = paste0(vcf.meta$SVTYPE, " chr", vcf.meta$CHROM, ":", vcf.meta$POS, " (len:", vcf.meta$SVLEN, " maf:", vcf.meta$MAF,")")
res_final = vcf.meta[,c("ID","sv_info","closest_gene")] %>% right_join(res.all, by = c("ID"="sv_id")) %>% dplyr::rename("pval" = "nom_p") %>% arrange(fdr)
save(res_final, file = paste0(data.dir,"/",studyDataID,"/association_results.RData"))# Load genotypes and phenotypes
load(paste0(data.dir,"ROSMAP_SV_4_association.RData"))
load(paste0(data.dir,"/",studyDataID,"/association_results.RData"))
phenotypes = phenotypes %>%
filter(study==studyID)
res_final = res_final %>% arrange(pval) %>%
filter(ID %in% valid_MAF_SVs) %>%
filter(ID %in% valid_HWE_SVs) %>%
filter(pheno %in% names(pheno_list))
res_final = add_LD_info(res_final)
data.table::fwrite(res_final, file = paste0(data.dir,"/",studyDataID,"/association_results.tsv.gz"))
createDT(res_final %>% head(1000))ci = 0.95
df_m = res_final %>%
group_by(pheno) %>%
mutate(observed = -log10(sort(pval)),
expected = -log10(ppoints(n())),
clower = -log10(qbeta(p = (1 - ci) / 2, shape1 = 1:n(), shape2 = n():1)),
cupper = -log10(qbeta(p = (1 + ci) / 2, shape1 = 1:n(), shape2 = n():1)),
lambda = median(qchisq(1 - pval, 1)) / qchisq(0.5, 1)) %>%
reframe(lambda = unique(lambda), lambda_label = sprintf("%s (λ = %.2f)", pheno, lambda)) %>% distinct()
data.table::fwrite(df_m, file = paste0(data.dir,"/",studyDataID,"/lambdas.tsv.gz"))
createDT(df_m)## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.5 data.table_1.15.4 vcfR_1.15.0 lubridate_1.9.3
## [5] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
## [9] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
## [13] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.13 ape_5.8 lattice_0.20-45 digest_0.6.36
## [5] utf8_1.2.4 R6_2.5.1 evaluate_0.24.0 highr_0.11
## [9] pillar_1.9.0 rlang_1.1.4 rstudioapi_0.16.0 vegan_2.6-6.1
## [13] jquerylib_0.1.4 Matrix_1.6-5 DT_0.33 rmarkdown_2.27
## [17] labeling_0.4.3 splines_4.1.2 pinfsc50_1.3.0 htmlwidgets_1.6.4
## [21] munsell_0.5.1 compiler_4.1.2 xfun_0.46 pkgconfig_2.0.3
## [25] mgcv_1.9-1 htmltools_0.5.8.1 tidyselect_1.2.1 fansi_1.0.6
## [29] permute_0.9-7 viridisLite_0.4.2 tzdb_0.4.0 withr_3.0.0
## [33] MASS_7.3-60.0.1 grid_4.1.2 nlme_3.1-153 jsonlite_1.8.8
## [37] gtable_0.3.5 lifecycle_1.0.4 magrittr_2.0.3 scales_1.3.0
## [41] cli_3.6.3 stringi_1.8.4 cachem_1.1.0 farver_2.1.2
## [45] bslib_0.8.0 generics_0.1.3 vctrs_0.6.5 ggeasy_0.1.4
## [49] RColorBrewer_1.1-3 tools_4.1.2 glue_1.7.0 hms_1.1.3
## [53] crosstalk_1.2.1 parallel_4.1.2 fastmap_1.2.0 yaml_2.3.10
## [57] timechange_0.3.0 colorspace_2.1-1 cluster_2.1.2 knitr_1.48
## [61] sass_0.4.9